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Abstract 

We propose a new method for the extraction of Rotation Measure from spectral 

polarization data. The method is based on maximum likelihood analysis and takes 

into account the circular nature of the polarization data. The method is unbiased 

and statistically more efficient than the standard \ 2 procedure. We also find that the 

method is computationally much faster than the standard \ 2 procedure if the number 

I/"") . of data points are very large. 

(N 

o 

1 Introduction 

o 

Polarizations of radio waves from cosmologically distant sources undergo Faraday rotation 
upon propagation through galactic magnetic fields. This effect provides very useful infor- 
mation about the magnetic fields in our galaxy as well as in the host galaxy (Vallee 1997; 



Or 



o 



Zeldovich 83). The amount of rotation is proportional to the magnetic field component par- 
allel to the direction of propagation of the wave and to the square of the wavelength A. The 
observed orientation of the linearly polarized component of the electromagnetic wave 0(X 2 ) 
can therefore be written as, 

e(\ 2 ) = e + (rm) x 2 

where the slope, called Rotation Measure {RM), depends linearly on the parallel compo- 
nent of the magnetic field and 6q is the intercept, also called the intrinsic position angle of 
polarization, IP A. 

RM and IP A can be determined by the standard \ 2 procedure. The extraction of RM is 
ambiguous since the observed polarization is defined only upto additions of nir, where n is an 
integer (see for example, Simard-Normandin et al. 1981; Zeldovich et al. 1983; Broten et al. 
1988). Therefore one finds many fits depending on the choice of nir for the polarizations at 
different wavelengths and it is necessary to find the RM which gives the absolute minimum 
value of x 2 ■ However, for small data sets choosing absolute minima is not reliable since it is 
often possible to make \ 2 arbitrarily small by making RM sufficiently large. It is therefore 
necessary to limit the range of RM while searching for the minimum \ 2 . 
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Another problem that also arises in many cases is that the observed polarization does 
not actually show a linear dependence on A 2 (Vallee 1980). This can arise if different regions 
of the source galaxy dominate at different wavelengths. It has, therefore been recommended 
that one should concentrate on a narrow region of the wavelength where the straight line 
may be a good fit to data. A methodology for the selection of limited wavelength region 
based on a unique physical regime has been suggested by Vallee (1980) and used by Broten, 
MacLeod & Vallee (1988) to compile a catalogue of unambiguous RM. 

In the present paper we propose an alternative method for the extraction of RM and 
IP A from data. This procedure is based on maximum likelihood for distributions defined 
to be invariant under angular coordinate changes. The von Mises (vM) distribution serves 
as a prototype for the statistical fluctuations for circular data. It is given by, (Mardia, 1972; 
Batschelet, 1981; Fisher, 1993) 



cxp 



kco$2(9-9) 



MB) = l —r7-, (1) 



7r/ («) 

where k is a parameter which measures the concentration of the population and 9 is the 
mean angle. The factor 2 has been inserted since the polarization angle is ambiguous by 
rm and not 2mr (Born and Wolf, 1980; for a detailed discussion on this point see Ralston 
& Jain, 1999). The von Mises distribution for circular data is in many ways the analoque 
of normal distribution for linear data (Mardia, 1972; Batschelet, 1981; Fisher, 1993). For 
example, a two dimensional normal distribution g(x,y) restricted to the circle x 2 + y 2 = 1 is 
equivalent to the von Mises distribution. For large k the von Mises distribution is strongly 
peaked at = 9. Keeping only the leading order term in 9 — 9 we find that it reduces to 
the normal distribution in this limit. The von Mises distribution has found applications in 
many physical examples involving circular data. For example, it has been found (Kendall 
and Young, 1984; Jain and Ralston, 1999) that the variable (3 — 9 — ip, where 9 is the IPA 
and ijj is the orientation axis of the radio galaxy, is well described by this distribution. 

We can obtain the distribution of the polarization angle by making the reasonable as- 
sumption that the electric field components follow a normal distribution (Brosseau, 1998), 

P(Ei,R2) = ^^exp f- £ E\R ld E 3 

We work in the circular polarization basis and E 1 , E 2 here refer to the right and left circularly 
polarized components of the electric field. Rij are the parameters of the distribution. We 
make a change of variables such that Ej = aj exp(i0j). We are interested in the distribution of 
the linear polarization angle 9, which in the circular polarization basis is equal to (9i — 62) /2. 
After integrating over the rest of the variables we obtain the MacDonald-Bunimovitch (MB) 
distribution p M b(9) (Brosseau 1998, MacDonald 1949, Bunimovitch 1949), 



Pmb(9) 



7T(1 - ^ 2 ) 3 / 2 



//sm _ V+ t:V + (1 - /" 2 ) 1/2 



(2) 



where // = £ cos(2# — 29). Here £ 2 — 1 _ |i? 12 | 2 / R X1 R 22 (0 < £ < 1) and 9 are the parame- 
ters of the distribution which measure the concentration and mean value of the population 
respectively. 
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Figure 1: Comparison of the normal distribution (er = 0.8) with von Mises and the MB 
distribution, Eq. 2. The parameters of the von Mises and MB were determined by requiring 
that the least square difference between them and the normal distribution is minimum. 

The von Mises as well as the normal distribution reduce to the Pmb{@) in the limit when 
the population is strongly concentrated near the mean value, i.e. in the limit when the 
parameter £ — > 1. However if the distribution is not too narrow these distributions differ 
significantly as shown in fig. 1. Furthermore the von Mises distribution and Pmb{9) are 
invariant under angular coordinate translations which is not the case for normal distribution. 
The MB distribution turns out to be more sharply peaked in comparison to von Mises. It 
is interesting that the distribution of f3, the difference between the polarization angle and 
galaxy orientation angle, was empirically found to be better described by a more sharply 
peaked function in comparison to von Mises (Jain and Ralston, 1999). 

In our likelihood analysis we find it convenient to first determine the rotation measure 
and IPA using the von Mises distribution. This is done by maximizing the log likelihood 
function numerically. This extremization turns out to be relatively easy since it is possible 
to analytically eliminate one of the variables. It is then only necessary to perform a one 
dimensional extremization numerically. The rotation measure and IPA determined by this 
procedure then serves as an initial guess for the analysis using the more complex MB dis- 
tribution, Eq. 2, for which we have to maximize the log likelihood function over both the 
variables numerically. Since the polarization angle is expected to follow this distribution our 
procedure will yield an unbiased and statistically efficient measure (see for example, Freund, 
1992; Cramer, 1958) of the RM in comparison to the standard method of y 2 minimization. 
The x 2 procedure is based on the gaussian distribution which is a not invariant under angular 
coordinate changes and hence cannot provide an optimal description for angular data. We 
point out that the standard theorems, which assure that the x 2 estimate is unbiased and of 
minimal variance, are not applicable to circular data. 

Besides the fact that the maximum likelihood procedure is more efficient we find that 
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it is also computationally less intensive in the case that the number of data points are 
large. In the case of standard x 2 procedure a fit has to be made for all possible choices 
of the polarization angle 9 obtained after addition by ±nn (n = 0, ±1, ±2, • • •). With the 
increase in the number of data points in the chosen wavelength range the total number of 
combinations increase very rapidly requiring very long computation time to exhaust all the 
possible choices. However, large number of data points are preferable in order to reduce the 
ambiguities in the extraction of RM and IP A. We find that the computational time for 
the current procedure increases very slowly with the number of data points. This is because 
for the case of von Mises distribution it only requires a search over one parameter. By 
explicit calculations we find that the program converges to the absolute minima relatively 
fast without requiring an unduly small search interval. This solution is then used as an input 
for the more realistic MB distribution (Eq. 2), which also turns out to converge very fast to 
the true solution. Hence the procedure can be very useful in extracting unambiguous RM 
from data. 

We apply our procedure to the spectral data given in the catalogue of linear polarization 
by Tabara & Inoue (1980). The catalogue contains position angle of polarization as a function 
of A 2 for 1510 sources. However, for many of these sources there exist only one or two data 
points and hence the extraction of RM is not possible. For the remaining 704 sources RM 
is listed in the catalogue. This is obtained by limiting the range of RM search to less than 
200 rad/m 2 except on the galactic plane (cot |6| < 2). For cot |6| > 2 the range is extended 
to 100 cot \b\rad/m 2 . For three of these 704 sources, data is listed in the catalogue at only 
two wavelengths. We ignore these three sources for our analysis and apply our procedure for 
the remaining 701 sources. 



2 Likelihood Method for Extracting RM 

Our procedure is based on the maximum likelihood analysis using a distribution for po- 
larization angle obtained under the assumption that the electric field vector is normally 
distributed. We find it useful in actual computations to first use a simpler form for the 
distribution function which is invariant under angular coordinate transformations. We as- 
sume this to be the von Mises distribution given in equation (1). Since the mean angle 
9 = RM\ 2 + 6> , the distribution can be written as, 

= exp [K cos2( 9 - fi MA3-M 
■kI {k) 

The best fit parameters RM, 9q can be obtained by maximizing the log likelihood of the 
spectral polarization data 6>(A 2 ). Let 9 t and Aj denote the measured polarization angle and 
the corresponding wavelength respectively. The log likelihood of the data is then given by, 

N N 

In £ = «• cos 2(8i - RMX 2 - 9 ) - N In tt - ]T In J (^) (4) 

i=l i=l 

where N is the total number of data points for a given source. We can estimate the parameter 
Ki from the error in the measurement of the polarization angle. For a circular variable 29i 
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the measure of concentration (Batschelet, 1981) is given by the length of the mean vector 



1 \ M 

-J £ cos 2 (fly-fc) (5) 



where refers to a particular measurement of Oi and M are the total number of polarization 
measurements taken at a particular value of Aj. For von Mises distribution 

where is the modified Bessel function of order v. Assuming that the error in the measure- 
ment of 9i is small, we expect Ki to be large and ~ 1. Expanding the Bessel function for 
large value of Hi we find that 

1 



2(1 -n) 

If the assumption that the errors are not small is not valid then it is better to solve Eq. 6 
numerically. Since the error in the measurement of is A6*j we find that r\ = cos(2A#j). 
Therefore 



2 (1 - cos(2A0i)) 

Substituting this in the equation (4) and maximising the log likelihood with respect to RM 
and 6*o gives the estimate of these parameters. This turns out to be equivalent to minimizing 

2 ^ l-cos2[9, t -9 -(RM) Af] 

Xcir ^ 1-COS2A0; {) 

i=i £ 

which is analogous to the standard x 2 used in the case of linear variables. For a large Ki 
we find that xlir follows the x 2 distribution with iV degrees of freedom, the same as that 
followed by the standard linear x 2 - We find by explicit calculations that the xlir * s close to 
the usual x 2 f° r most of the sources. Minimising with respect to 9 gives 

E i sing < /(l-cos2AgQ 
ta °" EiCosS i /(l-cos2A0 i ) 

with Bi — 2[6 i — (RM) Af]. Substituting this into equation (7) we can find the minimum of 
xtir by searching numerically over a single parameter RM. 

For the case of more realistic MB distribution, Eq. 2, we again maximize the log likeli- 
hood, 

N 

ln£ = ]Tm[p MB (0,)] 
i=i 

as a function of the parameters RM and 8$ in exact analogy as was done for the von Mises 
distribution. The main difference in the present case is that it is not possible to eliminate 
9 analytically. Hence this distribution requires a two dimensional extremization. We use 
the results obtained with von Mises distribution as an initial guess for the present case. The 
two dimensional extremization uses the program 'powell' (Press et al 1986) which gives very 
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fast convergence. We point out that the parameter £ in the present case can be estimated 
for each polarization measurement from the corresponding error in that measurement. Since 
the errors are quoted assuming gaussian statistics, we find the parameter £ such that the 
corresponding MB distribution is closest to the normal distribution for the measured polar- 
ization angle and its error. By being closest we mean that the square difference between 
the two distributions is minimum. By explicit calculation we find that for small values of 
the standard deviation parameter a used in normal distribution, 1 — £ 2 1.3<j 2 . In order 
to obtain an estimate of the goodness of fit by MB distribution we define the log likelihood 
difference 

A(ln£) = £ [\n\p MB (0)} - ln[p MB (^)]] . (8) 
i=i 

By making a Taylor expansion around 9{ = and keeping terms upto Of this reduces to the 
standard x 2 upto an overall constant equal to 6/1.3. The 1.3 in the denominator arises due 
to factor of 1.3 between a corresponding to gaussian statistics and £ which appears in MB 
distribution. We find by explicit calculations that A(ln£) is proportional to the standard 
X 2 only for very small errors in the polarization angle. For large value of x 2 it does not 
have a linear dependence on x 2 ■ lu giving results for the MB distribution we quote this 
log likelihood difference. The x 2 /dof = 1 for gaussian statistics, which represents the cutoff 
point for what is considered as a good fit, correspond in the present case to A(\n C)/ do f ~ 2. 
For x 2 '/dof > 1 we find that the A(\n C)/ do f increases very slowly with increase in x 2 ■ 

3 Results and Discussion 

In Fig. 2 we show Xcir as a function of RM for several representative cases. We find 
that xlir nas several local minima. These minima corresponds to the standard ambiguity 
in the determination of RM due to rm ambiguity in the polarization angle. For large 
number (> 10) of data points, we find that unique minimum is selected within the range 
abs(RM) < 1000 radians /m 2 which turns out to be always in agreement with the value listed 
in the Tabara and Inoue (1980). We see from Fig. 2 that the search interval in RM does not 
have to be unduly small in order that the program converge to the true minimum and hence 
does not require very long computational time. The procedure can be further optimized by 
using several well known routines, such as simulated annealing, genetic algorithms designed 
to find a global minimum efficiently 

As we mentioned earlier Tabara and Inoue (1980) resolve the n7r ambiquity by searching 
the minimum in x 2 over a limited range of RM which is selected on the basis of the position 
of the source with respect to the galactic plane. We also use the same range of RM in order 
to search for the global minima. In most cases we find that our results are in good agreement 
with the ones given in Tabara and Inoue (1980), which have been obtained using the standard 
linear fit. A randomly selected sample of results using both the distributions (von Mises and 
the MB) are shown in Table 1. The complete set of results including the errors in RM and 
IP A can be found at the website http : / / home.iitk.ac.in/~pkjain/ RM _DAT A. 

In many cases, however, we find that the global minima of x 2 C iri within the limited region 
of RM, is not the same as the one obtained using standard x 2 ■ The global extrema for the 
log likelihood of MB distribution also differs in some cases from the extrema of both the 
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Figure 2: A representative sample of results showing xti r as a function of RM. 
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xtir and X 2 ■ We find that for 16 sources the solution selected by xlir * s different from that 
selected by x 2 ■ For all except two of these cases the MB distribution prefers the solution 
obtained by using xlir- Overall the MB distribution gives results different from xlir and 
X 2 for 30 and 44 sources respectively. A sample of results illustrating the different minima 
selected by these different methods are given in Table 2. These differences arise since the 
different procedures prefer different nir combinations. Hence in this sense the determination 
of rotation measure by the standard x 2 method does have some bias, which is related to the 
well known nn ambiquity in the polarization angle. 

In Fig. 3 we show RM L — RM M b as a function of RM L where RM L is the rotation 
measure calculated by minimization of x 2 , i-e. the standard method, and RMmb is the result 
obtained by maximizing the log likelihood for the MB distribution. In figure 3 we choose the 
solution RM MB , among the different possible extrema of the log likelihood, such that RM MB 
is closest to RM L . In Fig. 4, we show the comparison between the two determinations by 
choosing the absolute minima for the log likelihood within the RM search interval. We 
also calculate the correlation between RMl — RMmb and RMl in order to investigate the 
existence of systematic difference between the two determinations. We find that if choose the 
log likelihood extrema such that it gives RM MB closest to RM L then there does not exist any 
correlation between RM L — RM MB and RM L . We find in this case that np 2 = 0.05, where 
p is the statistical correlation parameter. However if we choose the absolute minima for 
determination of RMmb we find that there does exist correlation between RMl — RMmb 
and RMl (np 2 = 26). Hence in this sense these two determinations of rotation measure 
show systematic difference, reflecting the existence of bias due to nn ambiquity. Similar 
calculation has been done for IP A also. The results in this case for the correlation between 
the vectors [sm(IPA L ), cos(IPA L )] and [sm(IPA L - IPA MB ),cos(IPA L - IPA MB )] are 
shown in tables 3 and 4 for the two choices of log likelihood extremum discussed above. It is 
clear from these tables that the IP A determined by minimization of x 2 does show systematic 
difference from that determined by maximizing log likelihood of the MB distribution. 

It has been found that there are few cases shown in Fig. 5 where the xlir P er degree 
of freedom is very large (> 30). We have verified that the corresponding fit using MB 
distribution is also quite close to the fit using xlir- m the graphs we prefer to show the 
results obtained using vM distribution because the corresponding statistic xlir has a more 
direct interpretation than the likelihood difference A(ln£) in the case of MB distribution. 
For the cases displayed in figure our determination of RM is quite close to Tabara and Inoue 
(1980) values. We have verified that the standard linear x 2 is also very large for these cases, 
as expected. Hence it seems that the straight line is not a good fit for these sources. The 
dominant reason for this behavior appears to be that different regions in the source, which 
in general have different polarization, may dominate at different wavelengths. In order to 
make an unbiased extraction of RM for these sources it may be necessary to use a limited 
range of wavelength (Vallee 1980) chosen on the basis of some physical characteristics of the 
source. Alternatively one can introduce higher order terms in A 2 in the fit. The procedure 
proposed in this paper can be easily extended to include such terms. 
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Table 1: A randomly selected sample of the results using the likelihood analysis (Circular 
fit) with the MacDonald-Bunimovitch (MB) and von Mises (vM) distributions and its com- 
parison with the standard linear fit. The RM values are given in units of rad/m 2 and the 
IP A are in degrees. The A(ln £)/ do f , xlir/dof and x 2 /dof are the measures of goodness 
of fit for the three different cases, dof stands for degrees of freedom. 
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Table 2: A sample of sources for which the global minima of x1i r ano ^ X 2 occur at different 
RM values. A single value for the circular fit for the source 1957+405 indicates that only 
one extrema was found. 
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Table 3: Correlation matrix np 2 AB between [sm(IPA L ), cos(IPA L )] and [sin(IPA L — 
IPA MB ),cos(IPA L — IPAmb)] for the case when the extremum of log likelihood of MB 
distribution is chosen such that it gives rotation measure nearest to that obtained by mini- 
mizing x 2 ■ IPAl and IPAmb are the estimates of the intercept of the fit determined using 
the linear method and the likelihood analysis using MB distribution respectively. 
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Table 4: Correlation matrix np 2 A B between [sm(IPA L ), cos(IPA L )] and [sin(/PA L — 
IPAmb), c os(IPAl — IPAmb)} for the case when the absolute extremum of log likelihood 
of MB distribution is chosen within the search range of the rotation measure. IPA^ and 
IPAmb are the estimates of the intercept of the fit determined using the linear method and 
the likelihood analysis using MB distribution respectively. 
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Figure 3: Comparison of the RM determined by minimization of x 2 (RMl) with that de- 
termined by use of the MB distribution (RMmb)- The units of RM are radian/m 2 . Here 
the RM MB solution is selected such that it gives results closest to the linear determination, 
even if the absolute maxima in the log likelihood occurs at a different location. The range in 
y-axis has been limited to ±40 for the sake of clarity. There were three data points outside 
this range which are not shown in the graph. 
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Figure 4: Comparison of the RM determined by minimization of \ 2 {RMl) with that de- 
termined by use of the MB distribution (RM MB ). The units of RM are radian /m 2 . Here 
the RM MB solution is selected by the absolute maxima of the log likelihood. 

4 Conclusion 

We have proposed a new method for the computation of rotation measures from spectral 
polarization data. The method takes into account the circular nature of polarization ob- 
servable and is based on the maximum likelihood analysis. It employs the MB distribution 
of polarization angle which can be derived by assuming that the electric field vector follows 
a normal distribution. Since the measured polarization angles are expected to follow this 
distribution the method will give an efficient and unbiased estimate of RM and IP A. A sim- 
pler procedure for determination of rotation measure, based on the von Mises distribution, 
is also presented. Due to its simplicity and invariance under circular coordinate transfor- 
mation, the von Mises distribution turns out to very useful for the purpose of providing an 
initial ansatz which can then be used for the more accurate RM determination using the 
MB distribution. For many purpose this initial ansatz may itself be sufficiently accurate. 
The von Mises distribution leads us to define a new statistic, x\ ir , which turns out to be 
very useful for handling circular data. In particular it provides us with an easy method to 
search over the RM range to find the absolute minima without requiring us to look at all 
possible nir combinations which can be added to the polarization angle. For this reason it 
turns out to be computationally very fast compared to the standard \ 2 if t ne data set is 
very large (N > 10). Since large data sets are required in order to obtain an unambiquous 
value of RM, the present method can be very useful. We further find that due to the nir 
ambiquity the rotation measure determined by the standard x 2 method shows systematic 
difference from the procedure using the MB distribution. This difference dissapears if we 
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insist that the same nir combination is used for the two determinations. However we find 
that systematic difference in the determination of IP A remains even if we require that the 
two methods use the same nn combination. The procedure can also be easily extended to 
extract higher order terms, i.e. nonlinear dependence on wavelength squared, which often 
arise due to source morphology and can play an important role in characterizing the source. 
Hence our procedure may provide us with a useful tool for analysing large polarization data 
sets, if they become available in future, for the purpose of extracting rotation measures and 
other information about polarizations from radio sources. 
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b. Source Name :: 0624-058 
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d. Source Name :: 1547-795 
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e. Source Name :: 1641+399 
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f. Source Name :: 1648+050 

Circular fit Linear fit 



70 
50 
30 
10 
-10 
-30 



i — i — i — i — r 



RM = 8.6 
IPA = 17.9 



J L 



J I I I L 



200 400 600 800 1000 

k (cm ) 



Q 



2 

o 



o 



60 
40 
20 

-20 



i — i — i — i — i — r 



RM = 5.8 
IPA = 19.0 



J I I I I L 



_1_ 



200 400 600 800 1000 

2 2 
X (cm ) 



Figure 5: (a) to (f) :: Comparison of the fit found for the cases where the xlir/dof is 
very large (> 30). In all the graphs the units of RM and IPA are rad/m 2 and Degrees 
respectively. 
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